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ABSTRACT 

We perform controlled numerical experiments to asses the effect of baryon mass loss 
on the inner structure of large galactic dark matter haloes. This mass expulsion is in- 
tended to mimic both the supernovae and AGN feedbacks, as well as the evolution of 
stellar populations. This study is meant in particular for precursors of massive Early 
Type Galaxies, wherein strong AGN feedback (often dubbed "QSO mode" in galaxy 
formation models) has been proposed to remove on a short timcscalc, of the order of 
a few dynamical times, a substantial fraction of their baryons. In a previous paper we 
evaluated the observational consequences (size increase) of this process on the galactic 
structure (Ragone-Figueroa & Granato 2011). Here we focus on the distribution of 
dark matter in the galactic region. It is shown that the inner region of the DM halo 
expands and its density profile flattens by a sizeable amount, with little dependence 
on the expulsion timescale. We also evaluate the effect of the commonly made approx- 
imation of treating the baryonic component as a potential that changes in intensity 
without any variation in shape. This approximation leads to some underestimates of 
the halo expansion and its slope flattening. We conclude that cuspy density profiles in 
ETGs could be difficult to reconcile with an effective AGN (or stellar) feedback during 
the evolution of these systems. 

Key words: galaxies: formation - galaxies: evolution - galaxies: elliptical and lentic- 
ular, cD - galaxies: haloes - quasars: general - method: numerical 
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1 INTRODUCTION 

A long standing puzzle of post recombination cosmology 
based on Cold Dark Matter, independently of the presence 
of cosmological constant or curvature in the adopted cosmo- 
logical model, is the so called core-cusp problem (for a recent 
review see De Blok 2010). Since the 90's, a long series of N- 
body (gravity only) simulations, dating back to Dubinsky & 
Carlberg (1991), has produced Dark Matter Haloes (DMH) 
whose inner density profile is reasonably well described by 
a power law par" with a ~ — 1, i.e. a cusp. This is at odd 
with several observations, suggesting a flat density profiles 
of DM in the inner region of real galaxies, i.e. a core. The 
case for cored DM density profiles is rather strong in dwarf 
and disc dominated Low Surface Brightness (LSB) galax- 
ies, wherein the dynamics of visible matter tracks safely the 
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largely dominating gravitational field of DM. Also the dy- 
namics of normal spiral galaxies is best interpreted by means 
of cored DM mass models (e.g. Salucci & Frigerio-Martins 
2009 and references therein), while the presence of a core 
in DMH of large elliptical galaxies is difficult to assess and 
controversial (e.g. Memola, Salucci & Babic 2011; Sonnen- 
feld et al. 2011; Tortora et al. 2011), since their central region 
is even more gravitationally dominated by the baryon com- 
ponent, and because of the less straight-forward to interpret 
orbital structure. 

A widespread idea is that the solution of the core-cusp 
problem should be searched for in the physics of baryonic 
matter, that can affect to some extent the distribution of 
DM in the inner region of haloes, and is not included in the 
aforementioned computations. 

However, the first baryonic process expected to occur 
worsen the problem. It consists in some further contrac- 
tion of the central region of the DMH, induced by the 
much stronger concentration-collapse of dissipative baryonic 
matter (e.g. Blumenthal et al. 1986; Gnedin et al. 2004; 
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Abadi et al. 2010, Gnedin et al. 2011). The exact impor- 
tance of this process is somewhat debated. Different numer- 
ical works found different results, and significantly less con- 
traction than that predicted by the simple analytic estimate 
by Blumenthal et al. (1986). In any case, this mechanism 
produces some further steepening of the inner profile, with 
respect to the prediction of gravity only cosmological simu- 
lations. 

Nevertheless, several subsequent (and more complex) 
baryon processes occur after this primary condensation, and 
can act in the opposite direction. Indeed, they have been of- 
ten evaluated specifically to seek solutions to the core-cusp 
problem (e.g. Navarro, Eke & Frenk 1996; Gnedin & Zhao 
2002; Read & Gilmore 2005; Tonini, Lapi & Salucci 2006; 
Mashchenko, Couchman & Wadsley 2006; Mashchenko, 
Wadsley & Couchman 2008; Pasetto et al. 2010; Governato 
et al. 2010; De Souza et al. 2011; Ogiya & Mori 2011; Inoue 
& Saitoh 2011; Martizzi et al. 2011; Pontzen & Governato 
2012; Maccio et al. 2012). Several of these papers have ad- 
dressed the effect of baryon ejection, due to SNae feedback, 
on the central DM density profile of less massive systems, 
i.e. mainly dwarfs or at most disk dominated galaxies. The 
results on its importance have been sometimes contradic- 
tory (for a discussion of possible causes, see Read & Gilmore 
2005). 

Much less attention has been paid so far to these ef- 
fects in larger systems, such as Early Type Galaxies (ETGs). 
The main reason is obviously that their DM content is less 
constrained by the data, as noted above (for a review see 
Buote & Humphrey 2012). Moreover, mechanisms related 
to SNae feedback are expected to be less effective in the 
deeper potential wells of large ETGs. As for the first point, 
it is worth noticing that recent observational evidence points 
either to a DM distribution shallower (Memola, Salucci & 
Babic 2011) or cuspier (Sonnenfeld et al. 2011; Tortora et 
al. 2012) than that predicted by ACDM gravity only sim- 
ulations. New techniques could soon provide constrains on 
the DM distribution in ellipticals (e.g. Pooley et al. 2011). 
As for the second point, in the past few years it has become 
common to consider in galaxy formation theory the feedback 
of AGN activity, which is likely to be much more effective 
in massive systems than that due to SNae (e.g. Silk & Rees 
1998; Fabian 1999; Granato et al. 2001, 2004; Benson et al. 
2003; Cattaneo et al. 2006; Monaco et al. 2007; Sijacki et al. 
2007; Somerville et al. 2008; Johansson et al. 2009; Ciotti, 
Ostriker & Proga 2009) . AGN feedback is widely considered 
the most promising mechanism for relaxing those tensions 
between galaxy formation models and observation, broadly 
related to the so called overcooling problem. This is the ten- 
dency of models to produce too massive and too blue (i.e. 
star forming) galaxies at low redshift, and to lock up in 
galaxies an excessive fraction of available baryons. On the 
contrary, many observations indicate that on average the 
more massive a galaxy, the earlier it stops its star forma- 
tion activity, a phenomenon referred to as downsizing (e.g. 
Brammer et al. 2011; for a critical discussion of the various 
manifestations of downsizing see Fontanot et al. 2009). De- 
spite the likely prominent role of AGN feedback in galaxy 
formation, so far its effect on the DM distribution has not 
been evaluated in detail (but see Peirani, Kay & Silk 2008 
and Duffy et al. 2010). 

One kind of AGN feedback (that sometimes called 



"QSO mode") is expected to eject on a small timescale (of 
the order of a few 10 Myr at most) the cold gas not yet 
converted into stars in star forming ETGs. In a recent work 
(Ragone-Figueroa & Granato 2011; henceforth Paper I) we 
have evaluated with aimed numerical simulations the im- 
portance of this process, concentrating our analysis on the 
baryon component of ETGs, in order to asses its possible 
contribution to the observed size evolution of ETGs (see e.g. 
Newmann et al. 2011 and references therein). In this paper 
we investigate instead in detail the effects on the profile of 
the DM halo. 

Since we are mainly interested in ETGs rather than less 
massive systems such as dwarfs, we will extend the stud- 
ies already published to regions of the parameter space not 
covered previously. Indeed, this is the first work where the 
initial conditions have been thought to get a configuration, 
after the loss of a substantial fraction of baryons previously 
condensed in the central region of the DM halo (i.e. the 
galactic region), consistent with our basic knowledge of the 
properties of local large ETGs (baryon to DM mass ratio, 
scalelengths, size as a function of stellar mass). By converse, 
in most studies the baryons totally disappear, or the remain- 
ing fraction is 5%, consistently with the very low baryon 
content of dwarf galaxies. This is not the case for ETGs, 
where the leftover baryons still dominate the potential wells 
in the central region of the DM halo. Also, in order to asses 
the maximal effect of baryon loss, in most studies the ini- 
tial mass ratio of condensed baryon to DM has been set 
close to the cosmic baryon fraction. However this appears 
too extreme, at least for ETGs. According to the cooling 
prescriptions adopted by semi analytic models and simula- 
tions of galaxy formation, no more than a few tens percent of 
the cosmic baryons had time to cool and condense in a large 
galactic DM halo, during the few Gyrs over which ETGs 
formed at z > 1.5 (e.g. Benson et al. 2001; Helly et al. 2003; 

Granato et al. 2004; Cattaneo et al. 2007; Viola et al. 2008). 
We point out explicitly that the results of previous numer- 
ical experiments, when featuring substantial differences in 
the initial and/or final mass ratios of baryon to DM, or in 
the assumed density profiles of the two components (often 
disks for the baryons), cannot be simply re-scaled to predict 
quantitatively our findings. 

Moreover, a technical difference with respect to most 
previous similar works is that often the baryon component 
has been treated as a fixed shape (rigid) potential. Here in 
general we do not adopt this approximation (since in Pa- 
per I we wanted actually to study the expansion of leftover 
baryons), but we evaluate its effect. 

The organization of the paper is straightforward: in Sec- 
tion [5] we describe the initial conditions and method for our 
simulations, whose result are presented in Section [3] The 
implications for our understanding of ETGs formation are 
discussed in Section 2] 



2 NUMERICAL METHOD AND SETUP 

The setup of the simulations used in this paper is similar 
to that described in Paper I, to which we refer for a few 
more details. The main difference is that here we run also 
simulations in which the initial conditions take into account 
the possible halo contraction due to baryon condensation. 
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Figure 1. The curve marked NFW TEO represents the typical density run of DM haloes produced by gravity only simulations (Navarro, 
Frcnk & White 1997), whilst ABA TEO and BLUM TEO are intended to keep into account the contraction resulting from galaxy 
formation, as described in Section [2] After setting these theoretical initial conditions (dashed lines in both panels) we let each system 
to evolve, for 0.4 Gyr, to an equilibrium configuration. These are considered the initial conditions (IC) for the numerical experiments 
before forcing baryon mass loss. The density and mass profiles arc shown in the left and right panel, respectively. In addition we show 
in the right panel the mass Hernquist profile corresponding to the baryonic component. 



Moreover, we evaluated the effect of the approximation of 
treating the baryons as a fixed shape potential. 

The purpose of the simulations is to investigate the 
evolution of collision-less particles (stars and DM) under a 
change of gravitational potential due to a loss of baryonic 
mass of the system. In general, the escaping mass can repre- 
sent either the gas which has not been converted into stars 
during the star forming phase of the spheroid (ejected by 
feedback driven galactic winds), or the mass lost from stars 
in form of stellar winds and SNae explosions (which is likely 
to escape the ETGs potential wells). In any case, we assume 
as given, and due to causes not included in our physical 
treatment (such as SNae and AGN feedbacks, or stellar evo- 
lution), the temporal dependence of this mass loss (Eq. llO[) . 
which we put by hand, and we simulate the ensuing dynam- 
ical evolution of collision-less mass distributions. Therefore 
we don't have to treat the gas dynamics. 

We used the public version of the code GADGET-2 
(Springel et al. 2005) to perform simulations with 10 6 and 
5 x 10 6 particles, in the gravity only mode. None of the 
presented results shows any noticeable difference in the two 
cases, which assures us that the mass resolution is sufficient 
for the purposes of the present study. Half of the particles 
are used to sample the baryonic and dark matter compo- 
nents respectively, with a softening of 0.007 and 0.35 kpc 
respectively. We checked that our results are not affected by 
significant variations in these choices. 

The density distribution of DM particles is initially (i.e. 
even before computing the effects of baryonic contraction; 
see below) assumed to follow the standard NFW (Navarro, 
Frenk & White 1997) shape 



PDM(r) 



Mvir.DM c 2 g(c) 



(1) 



47ri?3. r r(l + cr) 2 ' 
where M v i r ,DM is the halo virial mass in DM (the DM mass 



inside R v i v ), f — r/R v - 1T , c is the concentration parameter 
and g(c) = [ln(l + c) - c/(l + c)]" 1 . 

The virial radius R v [ r is by definition the radius within 
which the mean density is A v i T (z v i r ) times the mean matter 
density of the universe p u (z v ir) at virialization redshift z V i r : 



Rvii — 



1/3 



(2) 



The overdensity A v i r («), for a flat cosmology, can be ap 
proximated by 

(18tt 2 + 82a; - 39a; 2 ) 



Q(z) 



(3) 



where x = Q(z) — 1 and is the ratio of the mean matter 
density to the critical density at redshift z (Bryan & Norman 
1998). The corresponding mass distribution is written 



M D m(< r) = Af vir ,DM g(c) 



ln(l + cr) 



(1 + cf) 



(4) 



For the collisionless baryonic particles (representing the 
potential of both stars as well as gas before expulsion), we 
assume that, as a result of the assembly of the central galaxy, 
they settle on an Hernquist (1990) profile, which provides a 
reasonable description of stellar density in spheroids: 

Mb a 1 



2ir r (r + a) 3 



(5) 



where Mb is the total baryonic mass. The corresponding 
mass distribution is 



M B (<r) = M B (— r —\ 
\r + aj 



(6) 



so that the half-mass radius is related to the scale radius a 
by R1/2 — (1 + v2) a and, assuming a mass to light ratio 
independent of r, the effective radius is R e ~ 1.81a. 

As discussed in the introduction, the process of baryon 
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condensation in the centre of the DMHs is expected to con- 
tract to some extent the DM distribution. Since, at variance 
with Paper I, the target of this study is precisely to evaluate 
the effect of baryon expulsion on DMH density profiles, we 
performed also runs including an estimate of the contraction 
in the initial conditions, by adopting the Abadi et al. (2010) 
prescriptions. They found that a simple formula captures 
the average behavior of their simulations: 



Ms(t) 



M B ( t =Q) exp 0^ t 



(10) 



r / /r i = l+a[(M i /M / ) n -l] 



(7) 



where n and a are parameters (see below). This equation 
relates the initial (i.e. before baryon condensation) and final 
radii of the spheres containing the same amount of DM, to 
the total masses Mf and Mi within the same spheres. Since 

M£ = M f , DM (<r f ) + M fMr {<r f ) 

Mi M i:DM (< n) + Mi Mr (< n) 

where A// : da/(< r/) = Mi,Di\i(< n), by definition of r; 
and r/. Once it is assumed an initial density distribution for 
DM and baryons and a final density distribution for baryons, 
Eq. [7] is an implicit equation for rj given an r;, that can be 
solved numerically to obtain the final mass distribution of 
DM, M/,dm(< r). As stated by Abadi et al. (2010), their 
results are well described setting n = 2 and a — 0.3. More- 
over the contraction predicted by Blumenthal et al. (1986), 
as well as that found in simulations by Gnedin et al. (2004), 
are well described setting n = 1, and a = 1 or a = 0.73 
respectively. When exploiting Eq. we describe the initial 
mass distribution of DM and the final mass distribution of 
baryons with Eq. 2] and Eq. [B] Also, we assume that before 
contraction baryons follow the same density distribution of 
DM (Eq. [1]), rescaled by the cosmic baryon fraction we as- 
sume f b = Q. bar /n m = 0.17 (consistent with CMB WMAP-7 
constraints, Komatsu et al. 2009). 

In the following, unless otherwise specified, by dynam- 
ical time tdyn we mean the initial (i.e. before any mass loss 
and expansion) dynamical time, computed at baryon Ri/2- 



tdyn — 



/2 



2G (M b /2 + M dm (< R 1/2 )) 



1/2 



(9) 



Even though in this paper our focus is on the evolution of 
the DM component of the system, we maintain the same 
definition of dynamical time as in Paper I, since the effects 
under study are confined to the inner region of the DM dis- 
tribution. 

Given the density profiles, we obtain the ID velocity dis- 
persion by integrating the Jeans equation under the assump- 
tion of isotropic conditions. To generate initial conditions, 
we randomly populate the system with baryonic and DM 
particles, according to the density distributions Eq. [5] and 
Eq. [1] The particles velocities are randomly generated as- 
suming local Maxwellian distribution, with ID velocity dis- 
persion obtained by the solution of Jeans equation. Then, by 
numerically evolving the system for several dynamical times, 
we obtain a (quasi-)static statistical equilibrium. This latter 
step does not produce significant variations in the density 
profiles (see discussion of Fig. [TJ in Section [3] for more de- 
tails). 

Starting from this initial setup, we introduce a mass 
loss, by removing exponentially during an ejection time in- 
terval At a fraction 1 — e of the baryonic mass : 



The mass loss is practically attained by decreasing cor- 
respondingly in time the mass of the baryonic particles sam- 
pling the density field. After the end of the mass loss period 
At, we let the system to evolve till it reaches a new equilib- 
rium configuration. 

We wanted also to study the effect of a common ap- 
proximation done when evaluating the evolution of the DM 
distribution under a baryon mass loss, namely to treat the 
baryonic component as a potential that changes in inten- 
sity without any variation in shape (a rigid potential; e.g. 
Navarro et al. 1996; Ogiya & Mori 2011). For this purpose, 
we run also simulations in which the positions of baryonic 
particles were not updated in time. 

The reference value for the initial (i.e. before any mass 
loss) ratio of virial mass (total mass within the virial radius) 
to baryonic mass is M v i r /M B (t=o) = 20. This value is in 
keeping with the fraction 20%-40% of cosmic baryons that 
can cool and condense in the central region of a large galactic 
DM halo at z ^ 1.5, according to the prescriptions adopted 
by semi analytic models and simulations of galaxy formation 
(e.g. Benson et al. 2001; Helly et al. 2003; Granato et al. 
2004; Cattaneo et al. 2007; Viola et al. 2008). Moreover, 
assuming that after the initial condensation the halo looses, 
due to feedback induced galactic winds, a fraction between 
20% to 80% of this "galactic" baryonic mass, it is left with 
a baryon to DM content broadly consistent with estimates 
in the local universe for large galactic haloes (Monster et al. 
2010). 

We set M v i r = 1O 13 M0 in all simulations. Nevertheless, 
our results apply to different values of M v i r , provided that 
the ratios of scale radii and masses in the two components 
(DM and baryons) are not changed, and that the time is 
measured in units of dynamical time tdyn oc p^ 1 ^ 2 . 

We adopt a concentration parameter c = 4, a typi- 
cal value at galactic halo formation (see Zhao et al. 2003; 
Klypin, Trujillo- Gomez, & Primack 2010), and R vil ~ 170 
kpc (from Eq. [2] and Eq. [3] with M vir = 10 13 M Q and 

Zvir — 3) . 

Finally, we set a = 1.5 kpc (7i e ~ 2.7 kpc). We refer 
the reader to Section 3 of Paper I for motivations for this 
values. The initial (i.e. before mass loss and expansion, Eq. 
O dynamical time as defined above is t^yn ~ 5 Myr. 

In summary, the parameters affecting the results of our 
simulations are the ratio of mass between the total and 
baryonic components M v i r /M B ( t =o); the corresponding ra- 
tio of scale- lengths i? v ir/a; the fraction of baryonic mass 
lost (1 — e), and the time At over which the loss occurs. We 
performed simulations covering broad ranges of the latter 
two quantities, while in most runs we kept the former two 
at the fiducial values reported above. We checked however 
that none of our qualitative conclusion is affected by factor 
~ 2 variations of them (Paper I and some more discussion 
in Section [3}. 



3 RESULTS 

Fig.[T]shows the initial equilibrium density and mass profiles 
for the DM component, i.e. those adopted before forcing any 
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Figure 2. Comparison between the initial and final mass profiles. The latter is the new equilibrium configuration reached ~ 0.2 Gyr after 
the baryonic mass loss (see Fig. [4j . Each panel refers to simulations performed adopting a different initial density profile, as indicated by 
the inner label, which is plotted with the same color as in previous figure. ABA HALFSIZE IC (bottom right panel) is obtained as ABA 
IC, but halving the Hernquist scale radius of the baryonic component. The final profiles for e = 0.2 and e = 0.8 are in magenta and cyan 
respectively; solid lines are for At = Myr (instantaneous baryon expulsion) and dashed lines for At = 80 Myr (barely distinguishable 
for the lower mass loss e = 0.8). In all panels we plot for reference also the standard NFW profile (black solid line). 



loss of baryonic mass (solid lines) . The right panel also dis- 
plays the adopted standard Hernquist mass profile for the 
baryons. Here and in the following, the lower limit of the 
density and mass profile plots equals the adopted gravita- 
tional softening for the DM component, the density and the 
accumulated mass are evaluated using radial bins contain- 
ing ~ 500 DM particles. These initial profiles are obtained 
letting to evolve for several tens of dynamical times the com- 
posite system, where for the DM component we adopted the 
pure NFW mass profile (black dashed line), or the profiles 
implicitly given by Eq. [7J with values of the n and a param- 
eters adequate to reproduce the contraction found by Abadi 
ct al. (2010) in simulations (red dashed line) or that ana- 
lytically predicted by Blumenthal et al. (1986) (blue dashed 
line), as detailed in Section [5] We refer to these three cases 
as NFW, ABA and BLUM respectively. The first and the 
last should be regarded as extreme cases where the effect 



of contraction is negligible and maximal, respectively, while 
the intermediate case, ABA, is likely to be more realistic. 
The numerical evolution does not affect much the initial 
configuration of DM. For instance, the local power law in- 
dex of the density distribution increases in the inner ~ 10 
kpc by less than 0.1, making it slightly steeper. The effect 
decreases with increasing contraction of the DM profiles, 
i.e. from NFW to BLUM. Kazantzidis, Magorrian & Moore 
(2004) suggested that the local Maxwellian assumption em- 
ployed here (and in many other works) to set up the initial 
conditions, can lead to spurious results in the study of long 
term evolution of DM haloes. In particular they found that 
the central cusp is quickly significantly reduced. In our case 
we find a much less important, and opposite, trend. We at- 
tribute this to the presence of the dominating baryonic com- 
ponent in the center of the halo, not included in the study 
by Kazantzidis et al. Indeed, on one hand we run a pure 
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Figure 3. The same as the previous figure but comparing the initial and final density profiles. The different panels correspond to different 
ICs as labeled in the plot. 



DM test case, which confirmed their result on the reduction 
of the central cusp. On the other hand, we note that also 
Read & Gilmore (2005), who considered composite baryon 
plus DM systems as we did, found excellent agreement be- 
tween initial conditions generated with the Maxwellian ap- 
proximation and with the alternative procedure proposed by 
Kazantzidis et al. (2004). 

In addition to the three cases NFW, ABA and BLUM, 
and in order to study the dependence of our results on the 
parameters of the initial baryonic configuration, we consid- 
ered ICs in which we halved the scale radius of the baryonic 
component (henceforth ABA HALFSIZE) or we doubled its 
mass. We then correspondingly contract the DM halo as 
we did for ABA. The direction of these variations has been 
elected in order to increase the effect of mass loss with re- 
spect to the standard case. These two initial conditions are 
not shown in Fig. [T] since they are almost indistinguishable 
from ABA IC. 

Fig. [2] and FigO display the effect of baryon mass loss 
on the initial mass distributions described above. The four 
panels in each figure correspond to the ICs dubbed NFW, 



ABA, BLUM and ABA HALFSIZE. The new equilibrium 
configurations, shown in the figures, are reached typically a 
few tens of dynamical times (Eq. [9} after the end of the mass 
loss period (At), as can be appreciated in Fig. [4] Whenever 
the fraction of baryon lost is important (say £ 50%), the 
final DM profile is significantly flatter in the center than the 
initial one. If this is the case, the ensuing expansion more 
than counteracts the opposite contraction caused by baryon 
condensation, at least that found in recent numerical sim- 
ulations (Abadi et al. 2010), see right panels. The effect is 
somewhat weaker when the DM profile is initially more con- 
centrated, as expected since in this case the contribution of 
DM to the gravitational field in the center becomes more im- 
portant. For the same reason, the flattening effect in ABA 
HALFSIZE is slightly enhanced with respect to ABA. In- 
deed, in the former ICs, the contribution to the initial equi- 
librium of the baryons that we let to escape during mass loss 
was higher. To a lesser extent, the same happens doubling 
the initial baryon mass, but the difference is barely visible, 
so that we don't plot this case. 

The dependence on the timescale of the mass loss is 
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Figure 4. Time evolution, after the baryon mass loss, of the density (left panel) and mass (right panel) DM profiles for one of our 
simulations, namely the case ABA with instantaneous mass loss At = Gyr and e = 0.2. The new equilibrium is reached within ~ 0.2 
Gyr, i.e. within 40 dynamical times of the central region, Eq.[9] This holds true for all the other cases. 
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Figure 5. Comparison between the new equilibrium DM density (left panel) and mass (right panel) profiles after baryon mass loss, for 
e = 0.2, obtained letting the baryon particles to adjust their position according to the change of potential (magenta lines), or keeping their 
position fixed (green lines). The former is the treatment adopted in this work, while the latter corresponds to the common approximation 
of treating the baryons as a rigid potential. 



weak. This finding is at variance with respect to the result 
by Ogiya & Mori (2011), who found that if the loss is slow 
enough, the DM profiles are much less affected, up to the 
point that after a while it returns very close to the initial 
conditions. We have checked that this result is not due to 
the approximation used by Ogiya & Mori of treating the 
baryons as a fixed form potential (see below). This differ- 
ence could arise in part, but not entirely, from the more com- 
pact baryon distribution adopted by these authors. Actually, 
we find a somewhat more evident dependence on timescale 
when the baryon component is initially more concentrated 
(ABA HALFSIZE; lower right panels). It is also interesting 
to note that the expansion of the baryonic mass distribution 



is much more affected than that of the DM by the timescale 
of expulsion (see Paper I for details). 

In Fig. [4] we show the time evolution of the density 
(left panel) and the mass (right panel) profiles of a repre- 
sentative model, after the baryon mass loss. The chosen run 
corresponds to the case ABA with instantaneous mass loss 
At — Gyr and e = 0.2. The figure shows that the profile 
at remains stable after ~ 0.2 Gyr. The new equilibrium is 
reached within 40 dynamical times of the central region (Eq. 
O. This holds true for all the other simulated cases. 

A commonly made approximation in evaluating the ef- 
fect of baryon expulsion (or condensation) on the DM halo is 
to treat the baryonic component as a potential that changes 
in intensity without any variation in shape. Fig.[5]illustrates 
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that this approximation leads to some underestimate of the 
slope flattening and of its (weak) dependence on the expul- 
sion timescale. This is due to the fact that the mass loss 
actually causes a significant broadening of the baryon con- 
centration in the centre, which changes the shape of the po- 
tential and is more important for shorter timescales (Paper 
I). Therefore, in the full computation, the DM distribution 
is flattened not only by the decrease of gravitational force 
of baryons, but also by the (timescale dependent) outside 
dragging, due to the expansion of the leftover baryonic mat- 
ter. 



4 DISCUSSION AND CONCLUSIONS 

A well know general prediction of cosmological, gravity only, 
simulations is that DM haloes should have cuspy density 
profiles, essentially independently of the mass scalqj. Obser- 
vations at small to medium galactic scales (dwarf galaxies, 
disc dominated LSB galaxies, as well as normal spirals) have 
demonstrated that this is not the case. At cluster scales, 
where DM is gravitationally subdominant in the central re- 
gion, the situation is instead far from clear, with several 
claims for cored (e.g. Ricthler et al. 2011) as well as for 
cuspy (e.g. Zitrin et al. 2011) density profiles. The former 
mismatch is widely ascribed to the back-reaction of baryons, 
whose evolution to form galaxies is driven by various non- 
gravitational processes, on DM particles. The intermediate 
regime of large ETGs, the subject of our work, has received 
so far much less attention, both from the observational as 
well as the theoretical point of view. As for the former, again 
there are interpretations favoring a cored (e.g. Memola et 
al. 2011) as well as a cuspy (Tortora et al. 2010, Sonnen- 
feld et al. 2011) density profile. Firm conclusions are how- 
ever severely plagued by the degeneracy between the Initial 
Mass Function (IMF) and the DM density profiles (Treu et 
al. 2010). 

In the present paper, we have elucidated that an impor- 
tant gas removal during the early evolution of ETGs, should 
leave as a byproduct sizeable signatures also on the inner 
profile of their DM haloes. A similar ejection is required by 
most galaxy formation models aiming to explain the basic 
properties of these systems, such as their chemical proper- 
ties, low baryon content or luminosity function (e.g. Benson 
et al. 2003; Granato et al. 2004; Pipino, Silk & Matteucci 
2009; Duffy et al. 2010), and it is commonly, but not always, 
ascribed to QSO mode AGN feed-back. The DM density pro- 
file ends up to be significantly less concentrated than NFW, 
unless the prior (opposite) contraction generated by baryon 
collapse and condensation has been very efficient, and prob- 
ably unrealistic, i.e. closer to that estimated on the basis of 
approximate analytical treatments (e.g. Blumenthal et al. 
1986), than to that found in most cosmological simulations 
(e.g. Abadi et al. 2010; Gnedin et al. 2011). 

Moreover, it has been pointed out that stellar feedback 
can weaken the DM cusp in dwarf and spiral galaxies not 
only by means of gas removal from the galaxy potential 
well, but also as a consequence of oscillations of the po- 
tential generated by bulk motions of gas within the galaxy 

1 note that this feature has been found recently also in self-similar 
analytic models for the halo collapse (e.g., Lapi & Cavaliere 2011). 



(Mashchenko et al. 2006, 2008; Potzen & Governato 2012; 
Maccio et al. 2012). As for the AGN feedback, similar fluc- 
tuations could be induced, for instance, by non isotropic 
gas removal. Although, for simplicity and lack of theoretical 
understanding, models of galaxy formation treat AGN feed- 
back by means of isotropic sub-grid prescriptions, the case 
for preferential directions for the effect of AGN activity on 
its environment is actually strong. Indeed, all our knowledge 
of the AGN phenomenon points to a non isotropic struc- 
ture, including accretion discs, jets, ionizations cones and 
dust tori. More specifically, recent observational evidence 
directly indicates non isotropic quasar driven gas removal 
(e.g. Cano-Diaz et al. 2011). To explore properly this effect 
would require more complex numerical experiments, possi- 
bly including a treatment of gas dynamics. The effort seems 
somewhat premature at present, given the lack of under- 
standing of how AGN gas removal works, and we leave it 
for future investigation^. In any case, it is worth noticing 
that this possible additional process could even enhance the 
flattening effect of AGN feedback on the DM distribution 
estimated here. 

Martizzi et al. (2011) suggested that several mecha- 
nisms contribute to the formation of the ~ 10 kpc core in the 
Brightest Cluster Galaxy of a simulated Virgo-like galaxy 
cluster (M V ir ~ 10 14 Mq), found when (and only when) sub- 
grid models for the growth of SMBHs and the ensuing feed- 
back are included in their hydro-simulation. Our idealized 
numerical experiments elucidates and quantifies the impor- 
tant contribution of one of these mechanismsa, namely the 
ejection of baryonic matter (representing gas). 

In conclusion, cuspy density profiles in ETGs, tenta- 
tively inferred from some recent observation, could be diffi- 
cult to reconcile with an effective AGN (or stellar) feedback, 
in particular that believed to cause massive galactic winds 
during the early evolution of these systems, making them 
red and dead. 
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